Dans le cadre du module de statistiques appliquées, nous avions abordé les bases des statistiques et pris en main des outils permettant d’effectuer des opérations statistiques avancées grâce au langage R. Ainsi, notre choix s’est porté sur des données statistiques relatives au prestige (score attribuée) de différentes professions canadiennes.
Ce projet de statistiques appliquées consiste en l’utilisation d’outils R adaptés afin d’effectuer une modélisation statistique sur notre jeu de données, d’étudier la relation entre plusieurs variables en effectuant une régression linéaire simple ou multiple, ainsi que produire un rapport.
R est un langage de programmation dont le but est de pouvoir traiter et organiser des jeux de données afin de pouvoir y appliquer des tests statistiques plus ou moins complexes et se représenter ces données graphiquement à l’aide d’une grande variété de graphiques disponibles. RStudio est une application proposant un environnement de développement et des outils adaptés au langage et à l’environnement de programmation R.
La fonction principale de RStudio consiste à faciliter le développement d’applications en langage R. Pour ce faire, le programme dispose de nombreux outils qui vous permettent notamment de créer des scripts, compiler du code, créer des graphes, ainsi que de travailler avec divers jeux de données.
L’extension R markdown permet de générer des documents de manière dynamique en mélangeant texte mis en forme et résultats produits par du code R. Les documents générés peuvent être au format HTML, PDF, Word, et bien d’autres. C’est donc un outil très pratique pour l’exportation, la communication et la diffusion de résultats d’analyse.
Notre jeu de données comporte 102 individus décrits par 6 variables :
head(data) %>% paged_table()viz_edu <- data %>%
ggplot(aes(x = prestige, y = education, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(6, 16, 2)
)
viz_edu <- ggplotly(viz_edu)
viz_inc <- data %>%
ggplot(aes(x = prestige, y = income, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(0, 20000, 10000)
)
viz_inc <- ggplotly(viz_inc)
viz_women <- data %>%
ggplot(aes(x = prestige, y = women, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(0, 100, 25)
)
viz_women <- ggplotly(viz_women)
viz_census <- data %>%
ggplot(aes(x = prestige, y = census, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(2500, 7500, 2500)
)
viz_census <- ggplotly(viz_census)
subplot(style(viz_edu, showlegend = FALSE), style(viz_inc, showlegend = FALSE), style(viz_women, showlegend = FALSE), viz_census, nrows = 2, margin = 0.07)
Dans un premier temps, nous pouvons observer une forte relation linéaire entre prestige et les variables education et income, contrairement à women et census.
Dans un second temps, nous pouvons étudier d’éventuelles corrélations entre deux variables à l’aide d’un nuage de points. Ici, nous utilisons la fonction ggpairs afin de pouvoir en observer plusieurs en même temps.
generate_scatterM <- function(data, mapping){
data %>% ggplot(mapping = mapping) + geom_point(size = 0.8) +
geom_smooth(method = loess, color = "red", size = 0.85, se = FALSE) +
geom_smooth(method = lm, color = "blue", size = 0.85, se = FALSE)
}
scatterM <- data %>% ggpairs(columns = 1:5, lower = list(continuous = generate_scatterM))
scatterM <- ggplotly(scatterM)
scatterM
Avec les observations précédentes, nous pouvons déjà écarter certaines linéarités et choisir lesquelles sont intéressantes à étudier. Ici, nous reprenons donc la relation entre la variable prestige (il s’agit d’un score de prestige relatif à la profession) et la variable education (qui reflète le niveau d’étude). Cette fois nous utilisons ggplot2 pour une meilleure représentation.
education <- data %>%
ggplot(aes(x = education, y = prestige)) +
geom_point() +
geom_smooth(method = loess, se = T) +
geom_smooth(method = lm, color = "red", se = F) +
theme_bw() +
scale_x_continuous(
breaks = seq(6, 16, 2)
) +
scale_y_continuous(
breaks = seq(20, 80, 20)
)
scatter_edu <- ggplotly(education)
scatter_eduLa droite bleue est définie par la méthode des moindres carrés (MSE), il s’agit d’une droite de régression linéaire entre les variables education et prestige.
La courbe rouge indique la tendance globale entre ces deux variables, il s’agit d’une courbe de régression de type lowess. Les deux courbes extérieures représentent l’intervalle de confiance de cette courbe de régression.
On constate que que la droite de régression est presque toujours comprise dans l’intervalle de confiance, l’hypothèse de linéarite entre les variables education et prestige est donc acceptable.
Nous testons aussi la relation entre les variables income et prestige.
income <- data %>%
ggplot(aes(x = income, y = prestige)) +
geom_point() +
geom_smooth(method = loess, se = T) +
geom_smooth(method = lm, color = "red", se = F) +
theme_bw() +
scale_x_continuous(
breaks = seq(0, 25000, 5000)
) +
scale_y_continuous(
breaks = seq(20, 100, 20)
)
scatter_inc <- ggplotly(income)
scatter_inc
Ici, en regardant la forme du lien entre la variable entre les deux variables, on s’aperçoit que la droite de régression suis bien moins l’intervalle de confiance de la courbe lowess. l’hypothèse de linéarité est alors plus critiquable.
Pour déterminer la droite de régression, on ajuste un modèle linéaire simple aux données, à l’aide de la fonction « lm ».
lin_model_edu <- lm(prestige ~ education, data = data)
lin_model_edu##
## Call:
## lm(formula = prestige ~ education, data = data)
##
## Coefficients:
## (Intercept) education
## -10.732 5.361« Intercept » correspond ici à l’ordonnée à l’origine, le « b » de notre droite, et le « x » est la pente de la droite, ce qui correspond au « b » dans notre notation.
L’équation, de notre droite est donc \(y = -10.732 + 5.361x\).
Il existe différentes manières, entre diagrammes et tests statistiques, visant à évaluer le lien linéaire entre deux variables. Ce lien est dit « significatif » s’il remplit certaines conditions. En effet, les résidus doivent être indépendants, distribués selon une loin Normale de moyenne 0 et de façon homogène (variance constante).
En premier lieu, il faut évaluer l’hypothèse d’indépendance des résidus. À l’aide d’un lag plot, il possible de mettre en évidence la présence d’une auto-corrélation.
bacf <- acf(residuals(lin_model_edu), plot = F)
bacfdf <- with(bacf, data.frame(lag, acf))
conf.level <- 0.95
ciline <- qnorm((1 - conf.level)/2)/sqrt(length(residuals(lin_model_edu)))
lag_plot <- bacfdf %>%
ggplot(aes(x = lag, y = acf)) +
geom_hline(aes(yintercept = 0)) +
geom_segment(aes(xend = lag, yend = 0)) +
geom_hline(aes(yintercept = ciline), linetype = 3, color = 'darkblue') +
geom_hline(aes(yintercept = -ciline), linetype = 3, color = 'darkblue') +
theme_bw()
lag_plot <- ggplotly(lag_plot, tooltip = c("lag", "acf"))
lag_plot
Sur le graphique ci-dessus, chaque trait correspond à un lag, ou coefficient de corrélation entre les résidus de chaque point, les pointillés bleus quant à eux sont les intervalles de confiance du coefficient de corrélation égale à 0.
En s’appuyant sur ce graphique, nous pouvons observer qu’une auto-corrélation significative est présente pour les lags 1 à 3, 5, 7 et 12 par exemple.
Afin d’aller plus loin et de détecter une éventuelle auto-corrélation des erreurs d’ordre 1 (lag = 1), il est possible d’emloyer le test de Durbin-Watson.
L’intérêt d’avoir recours à ce test est de déterminer si le modèle est perfectible. Souvent, l’autocorrélation s’observe sur les résidus d’une modélisation de série chronologique. Toutefois, si le type de modèle est bien choisi (en l’occurrence une régrassion linéaire simple), il peut exister une véritable autocorrélation entre les observations. En revanche, si aucune autocorrélation n’a de raison d’être, il faut chercher d’autres variables candidates.
durbinWatsonTest(lin_model_edu)## lag Autocorrelation D-W Statistic p-value
## 1 0.2752512 1.43917 0.002
## Alternative hypothesis: rho != 0
Le test de Durbin-Watson indique qu’il existe une auto-corrélation significative entre les résidus d’une ligne du tableau de données et ceux de la ligne suivante.
La deuxième condition susmentionnée concerne la distribution des résidus selon une loi Normale de moyenne 0. La pertinence de l’ajustement d’une distribution donnée à un modèle peut être représentée graphiquement grâce à un diagramme quantile-quantile, ou QQ plot.
qq_edu <- lin_model_edu %>%
ggplot(aes(sample = rstandard(.))) +
stat_qq_line(color = "red", linetype = "dashed") +
stat_qq(size = 1.2) +
theme_bw() +
scale_x_continuous(
breaks = seq(-2, 2, 1)
) +
scale_y_continuous(
breaks = seq(-3, 3, 1)
) +
labs(title = "Normal Q-Q",
x = "Theoretical Quantiles\nlm(prestige ~ education)",
y = "Standardized residuals")
qq_edu <- ggplotly(qq_edu)
qq_edu
Nous constatons que les résidus sont relativement bien alignés le long de la droite figurant sur le graphique, nous pouvons donc en conclure que la distribution de cette série de données suit une loi normale.
En statistiques, le test de Shapiro–Wilk teste l’hypothèse nulle selon laquelle un échantillon de données est issu d’une population normalement distribuée. Ce test peut être employé afin d’évaluer l’hypothèse de normalité des résidus.
shapiro.test(residuals(lin_model_edu))##
## Shapiro-Wilk normality test
##
## data: residuals(lin_model_edu)
## W = 0.98065, p-value = 0.1406On considère que cette hypothèse est rejetée si la p-valeur est inférieure à \(0,05\), mais dans le cas présent le test accepte la normalité.
Il est essentiel d’appliquer des tests statistiques tels que celui de Shapiro-Wilk. Un exemple parlant serait celui d’une régression linéaire simple entres les variables income et prestige.
lin_model_inc <- lm(prestige~income, data = data)
lin_model_inc##
## Call:
## lm(formula = prestige ~ income, data = data)
##
## Coefficients:
## (Intercept) income
## 27.141176 0.002897qq_inc <- lin_model_inc %>%
ggplot(aes(sample = rstandard(.)))+
stat_qq_line(color = "red", linetype = "dashed") +
stat_qq(size = 1.2) +
theme_bw() +
scale_x_continuous(
breaks = seq(-2, 2, 1)
) +
scale_y_continuous(
breaks = seq(-3, 3, 1)
) +
labs(title = "Normal Q-Q",
x = "Theoretical Quantiles\nlm(prestige ~ income)",
y = "Standardized residuals")
qq_inc <- ggplotly(qq_inc)
qq_inc
L’alignement des résidus est légèrement moins précis que sur le QQ plot entre les variables education et prestige, mais il semble plutôt satisfaisant. On ne constaste aucune valeur aberrante.
shapiro.test(residuals(lin_model_inc))##
## Shapiro-Wilk normality test
##
## data: residuals(lin_model_inc)
## W = 0.97169, p-value = 0.02729Cependant, le test de Shapiro-Wilk renvoie une p-valeur inférieure à \(0,05\) et rejette donc l’hypothèse de normalité.
Le concept d’homoscédasticité est utilisé dans le contexte de la régression linéaire pour décrire le cas où la variance des erreurs du modèle est la même pour toutes les observations. Autrement dit, les variances sont homogènes et les erreurs identiquement distribuées. Il s’agit de l’une des propriétés fondamentales du modèle de régression linéaire.
Afin de visualiser cela, nous avons besoin des valeurs prédites par le modèle (fitted values).
Dans notre cas, le graphique représentera les valeurs de prestige prédite par le modèle pour les valeurs de education présentes dans les données.
sc_loc_edu <- lin_model_edu %>%
ggplot(aes(fitted(.), sqrt(abs(rstandard(.))))) +
geom_point(size = 1.2) +
geom_smooth(method = loess, se = FALSE) +
theme_bw() +
scale_x_continuous(
breaks = seq(30, 70, 10)
) +
scale_y_continuous(
limits = c(0, 1.5),
breaks = seq(0, 1.5, 0.5)
) +
labs(title = "Scale-location",
x = "Fitted values\nlm(prestige ~ education)",
y = "sqrt |Standardized residuals|")
sc_loc_edu <- ggplotly(sc_loc_edu)
sc_loc_edu
Le graphique nous montre que les résidus sont répartis de façon homogène le long du gradient des valeurs prédites. Ceci est mis en évidence par la courbe de régression locale qui est quasiment plate.
Le test de Breusch-Pagan permet de tester l’hypothèse d’homoscédasticité du terme d’erreur d’un modèle de régression linéaire. Si la variance est constante, alors on a de l’homoscédasticité ; en revanche, si elle va varie, on a de l’hétéroscédasticité.
ncvTest(lin_model_edu)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.6327545, Df = 1, p = 0.42635La p-valeur associée au test est supérieure au seuil de \(0.05\), ainsi l’hypothèse d’homoscédasticité est acceptée.
Afin de vérifier l’hypothèse selon laquelle la relation entre les variables education et prestige est linéaire, ainsi que l’homoscédasticité, il est possible de générer un residuals vs. fitted plot.
res_v_fit_edu <- lin_model_edu %>%
ggplot(aes(fitted(.), residuals(.))) +
geom_point() +
stat_smooth(method="loess", se = FALSE) +
geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
theme_bw() +
scale_x_continuous(
breaks = seq(30, 70, 10)
) +
scale_y_continuous(
limits = c(-30, 20),
breaks = seq(-30, 20, 10)
) +
labs(title = "Residuals vs. Fitted values",
x = "Fitted values\nlm(prestige ~ education)",
y = "Residuals")
res_v_fit_edu <- ggplotly(res_v_fit_edu)
res_v_fit_edu
Il n’y a pas de schéma clair dans les données, ni d’aberrations évidentes. Les résidus sont, à peu de choses près, uniformément distribués le long de la ligne 0. La droite de régression est bien adaptée aux données, par conséquent l’hypothèse de linéarité est acceptable.
En s’appuyant sur les différentes représentations graphiques ainsi que les tests effectués précédémment, on peut en conclure que les propriétés fondamentales de la régression linéaire (indépendance des résidus, normalité, homoscédasticité et linéarité) sont satisfaites.
Partant de ce postulat, les résultats de la régression linéaire entre les variables education et prestige peuvent être interprétés.
La fonction summary nous permet d’avoir un bon aperçu du modèle. Elle affche les coefficients estimés, leur écart-type, et la valeur de la statistique t de Student ainsi que la p-valeur (probabilité que le coefficient soit significativement différent de zéro) associées à chaque coefficient. Sont aussi présentés le R2 et R2 ajusté, ainsi que la statistique F de Fisher (testant la significativité globale des variables), son degré de liberté, et la p-valeur associée.
summary(lin_model_edu)##
## Call:
## lm(formula = prestige ~ education, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.0397 -6.5228 0.6611 6.7430 18.1636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.732 3.677 -2.919 0.00434 **
## education 5.361 0.332 16.148 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared: 0.7228, Adjusted R-squared: 0.72
## F-statistic: 260.8 on 1 and 100 DF, p-value: < 2.2e-16
La vue du modèle est divisée en 2 grandes parties : Residuals et Coefficients.
Dans la première partie, nous avons la norme des résidus. Nous remarquons que la médiane est proche de 0, et que les valeurs absolues des 1er et 3e quartiles (1Q et 3Q) sont très proches. Cela porte à croire que les résidus sont distribués selon une Normale, ce qui est le cas au vu des résultats précédemment obtenus.
La seconde partie Coefficients est un tableau à deux lignes. Nous nous intéressons d’abord aux deux dernières colonnes. La valeur absolue de t calculée est bien supérieure à la valeur théorique déterminée, alors on rejette H0. Par conséquent, le coefficient est significativement différent de zéro et il est interprétable.
Dans le cas d’un modèle linéaire, le signe du coefficient associé à une variable indique le sens de l’effet de cette variable sur celle à expliquer. Sur la première colonne, le coefficient de la pente pour education est de \(5.361\), par conséquent augmenter le niveau d’éducation aura tendance à augmenter le score de prestige. Pour une unité supplémentaire concernant le niveau d’éducation, le prestige augmente de \(5.361\) unités (= coefficient de la pente).
La qualité de la prédiction d’une régression linéaire est mesurée par le coefficient de détermination R2. Plus il est proche de 1, et plus le rapport de convenance entre le modèle et les données est fort. Toutefois, ce coefficient R2 coît significativement avec avec le nombre de variables explicatives, c’est ici que le R2 ajusté intervient. En effet, il tient compte du nombre de variables et fournira une valeur plus pertinente.
Notre modèle correspond à une régression linéaire simple, donc les coefficients R2 et R2 ajusté sont très proches. L’un comme l’autre constitue une valeur pertinente. En s’appuyant sur l’avant dernière ligne, nous constatons qu’il est approximativement de \(0.72\), ce qui atteste d’un modèle fiable.
De plus, sur la dernière ligne la p-valeur est très basse et inférieure à \(0.05\), ainsi le lien linéaire entre les variables education et prestige est significatif.
Résidus de la régression :
data$residuals <- residuals(lin_model_edu)
head(data) %>% kable()| education | income | women | prestige | census | type | residuals |
|---|---|---|---|---|---|---|
| 13.11 | 12351 | 11.16 | 68.8 | 1113 | prof | 9.250875 |
| 12.26 | 25879 | 4.02 | 69.1 | 1130 | prof | 14.107621 |
| 12.77 | 9271 | 15.70 | 63.4 | 1171 | prof | 5.673573 |
| 11.42 | 8865 | 9.11 | 56.8 | 1175 | prof | 6.310758 |
| 14.62 | 8403 | 11.68 | 73.5 | 2111 | prof | 5.855950 |
| 15.64 | 11030 | 5.13 | 77.6 | 2113 | prof | 4.487854 |
Prédictions du modèle de régression pour les valeurs observées de la variable prédictive :
data$fitted <- fitted(lin_model_edu)
head(data) %>% kable()| education | income | women | prestige | census | type | residuals | fitted |
|---|---|---|---|---|---|---|---|
| 13.11 | 12351 | 11.16 | 68.8 | 1113 | prof | 9.250875 | 59.54913 |
| 12.26 | 25879 | 4.02 | 69.1 | 1130 | prof | 14.107621 | 54.99238 |
| 12.77 | 9271 | 15.70 | 63.4 | 1171 | prof | 5.673573 | 57.72643 |
| 11.42 | 8865 | 9.11 | 56.8 | 1175 | prof | 6.310758 | 50.48924 |
| 14.62 | 8403 | 11.68 | 73.5 | 2111 | prof | 5.855950 | 67.64405 |
| 15.64 | 11030 | 5.13 | 77.6 | 2113 | prof | 4.487854 | 73.11215 |
Matrice de variance-covariance des paramètres du modèle :
vcov(lin_model_edu) %>% kable()| (Intercept) | education | |
|---|---|---|
| (Intercept) | 13.520978 | -1.1835055 |
| education | -1.183506 | 0.1102162 |
Nous pouvons obtenir des réponses prédites par le modèle de régression, pour des valeurs de la variable prédictive n’ayant pas encore été observées.
predictions <- tibble(education=c(9.21, 11.07, 14.62))
predict(lin_model_edu, newdata = predictions, interval = "confidence") %>% kable()| fit | lwr | upr |
|---|---|---|
| 38.64170 | 36.58966 | 40.69374 |
| 48.61293 | 46.81134 | 50.41452 |
| 67.64405 | 64.52387 | 70.76423 |
Ici, nous avons réalisé des prédictions pour différents niveaux d’éducation, respectivement \(9.21\), \(11.07\) et \(14.62\) choisis arbitrairement.
pred_interval <- predict(lin_model_edu, interval = "prediction")
data <- cbind(data, pred_interval)linear_regression <- data %>%
ggplot(aes(y = prestige, x = education)) +
geom_point(size = 1.2) +
geom_smooth(color = "red", method = "lm", fill = "red") +
geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
geom_line(aes(y = upr), color = "red", linetype = "dashed") +
annotate("text", x = 9, y = 80,
label = paste0("prestige = ", coef(lin_model_edu)["(Intercept)"] %>% round(3),
" + ", coef(lin_model_edu)["education"] %>% round(3), " * education")) +
theme_bw() +
scale_x_continuous(
breaks = seq(6, 16, 2)
) +
scale_y_continuous(
breaks = seq(25, 75, 25)
) +
labs(title = "Linear Regression",
x = "Education", y = "Prestige")
linear_regression <- ggplotly(linear_regression)
linear_regressionLe modèle de régression linéaire est aussi bien utilisé pour chercher à prédire un phénomène que pour chercher à l’expliquer. Après avoir estimé un modèle de régression linéaire, nous pouvons prédire quel serait le niveau de y pour des valeurs particulières de x. Il permet également d’estimer l’effet d’une ou plusieurs variables sur une autre en contrôlant par un ensemble de facteurs. En apprentissage statistique, la méthode de régression linéaire est considérée comme une méthode d’apprentissage supervisé utilisée pour prédire une variable quantitative. Dans cette perspective, nous entraînons généralement le modèle sur un échantillon d’apprentissage et testons ensuite les performances prédictives du modèle sur un échantillon de test.
Enfin, le projet a été réalisé dans le cadre du cours de statistiques appliquées. Celui-ci convoque des notions en statistiques abordées durant cette année et fait le lien avec notre module d’apprentissage supervisé. Afin de concevoir nos systèmes prévisionnels, comprendre et interpréter les effets ou relations entre plusieurs variables d’un jeu de données est essentiel.